1 Présentation du document
2 Objectifs
3 Import des données
4 Préparation des données
5 Étude auto-corrélation et colinéarité
- 5.1 Réalisation d’une matrice de corrélation
- 5.2 Calcul de la Variance Inflaction Factor (VIF)
6 Comparaison des valeurs des communes fusionnantes vis-à-vis des valeurs des communes inchangées
7 Régression logistique classique
8 Étude des résidus (y compris avec indice de Moran pour une auto-corrélation spatiale des résidus ?)
9 Sélection step by step, Forward stepwise ou backward stepwise
Bibliographie

1 Présentation du document

NB : Ce script est un document de travail.

Il étudie les communes nouvelles en proposant un modèle pour observer si certaines variables sont déterminantes dans le fait de créer une commune nouvelle ou non.

Il est mis à disposition dans une logique de science ouverte.

Ce travail s’inscrit dans le cadre d’une étude plus générale sur les communes nouvelles :

https://cv.hal.science/gabriel-bideau

Licence CC-BY-NC-SA.

Il est possible d’accéder au code de ce Markdown ici : https://gbideau.github.io/CN_regression/regressions.Rmd

Les données utilisées pour jouer le code sont regroupées ici : https://gbideau.github.io/CN_data/

Ne pas hésiter à contacter l’auteur (gabriel.bideau@gmail.com) pour toute question.

2 Objectifs

L’objectif de ce script est de tester différents modèles de régressions logistiques permettant de modéliser le phénomène des communes nouvelles.

Ce travail s’inscrit dans le cadre de la thèse de G. Bideau sur les communes nouvelles, sous la direction de R. Le Goix.

Plan envisagé au départ :

Import de différentes données (par blocs complets)
Prévoir des tests sur des sous-ensembles (Maine-et-Loire, Savoies, Ouest…)Il est important voire intéressant de cartographier ou de représenter graphiquement les résidus. Cela permet potentiellement de valider ou invalider le modèle mais également de réfléchir aux espaces où ce modèle fonctionne le mieux ou le moins bien.
Calculs colinéarité, corrélations
Méthode de régression classique (copier-coller analyse budgets)
Forward stepwise ou backward stepwise
Cartographie des résidus (y compris avec indice de Moran pour une auto-corrélation spatiale des résidus ?)

3 Import des données

3.1 Géométries

On n’importe au départ que les données à la géométrie 2011, mentionnant les futures communes nouvelles de rattachement.

# Import des géométries 2011 avec les fusions uniquement jusqu'au 1er janvier 2022
geom2011 <- st_read("Archives/data_prep 2011-2022(01)/data/geom.gpkg", layer = "geom2011", quiet = TRUE) 
geomfus2011 <- st_read("Archives/data_prep 2011-2022(01)/data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) 
dep <- st_read("Archives/data_prep 2011-2022(01)/data/geom.gpkg", layer = "dep", quiet = TRUE)
geom_new <- st_read("Archives/data_prep 2011-2022(01)/data/geom.gpkg", layer = "geom_new", quiet = TRUE) 
geomCN_new <- st_read("Archives/data_prep 2011-2022(01)/data/geom.gpkg", layer = "geomCN_new", quiet = TRUE)  

# Import des géométries 2011 avec les fusions uniquement jusqu'au 1er janvier le plus récent (2024 au moment de l'édition de ce script)
geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE) 
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) 
dep <- st_read("data/geom.gpkg", layer = "dep", quiet = TRUE)
geom_new <- st_read("data/geom.gpkg", layer = "geom_new", quiet = TRUE) 
geomCN_new <- st_read("data/geom.gpkg", layer = "geomCN_new", quiet = TRUE)

3.2 Base DAC

Base DAC [bideau2022b]

# Import des données de la base DAC
# Si on souhaite importer les données comportant toutes les fusions uniquement jusqu'au 1er janvier 2022 (pour avoir également les données budgétaires) :
load("Archives/data_prep 2011-2022(01)/data/refdata.Rdata")
# Mais comme la régression s'appuie uniquement sur les données 2011, on peut prender la base DAC la plus récente
load("data/refdata.Rdata")

variables_RT <- colnames(df2011[which(stringr::str_detect(colnames(df2011), "_RT") == TRUE)])

3.3 Budgets

Cf. article dans la RERU, à paraître en 2024.

# Import des données budgétaires pour l'année 2011
load("data/refdata_budgets_bloccommunal_tot_2011-2022.Rdata")
rm(df_budgets_new, budgets_com_2011, budgets_com_2022)

# Fusion des données budgétaires avec les autres données
# On liste les variables renseignées dans les différents df
colbudgets2011 <- colnames(df_budgets_2011)

coldf2011 <- colnames(df2011)

# On liste les variables présentes seulement dans budgets (celles qu'il va falloir y prélever)
colgarder2011 <- setdiff(colbudgets2011, coldf2011)

# On fusionne les données de data_prep et des budgets, en ne gardant pour ces dernières que les données précédemment identifiées
df2011 <- merge(df2011,
                     df_budgets_2011[, c("CODGEO", colgarder2011)],
                     by = "CODGEO", all.x = TRUE, all.y = TRUE)

# Toutes les variables en pourcentages qui sont des variables budgétaires
variables_budg_prct <- colnames(df2011[which(stringr::str_detect(colnames(df2011), "_prct") == TRUE)])
variables_budg_prct <- variables_budg_prct[variables_budg_prct!= "geometry"] # On supprime la variable "geometry" qui est intégrée sans que cela soit voulu


# Suppression des données inutiles (déjà importées ou seulement de transition)
rm(df_budgets_2011, colbudgets2011, colgarder2011)

3.4 Données électorales 2012

PR2012_T1 <- read.table("data/elections/PR2012_T1.csv", sep="\t", colClasses = c(rep("character", 2), rep("numeric", 40)), head = TRUE, stringsAsFactors = TRUE, dec =",")
colnames(PR2012_T1)[3:42] <- paste0("PR2012_T1_", colnames(PR2012_T1)[3:42])

# PR2012_T2 <- read.table("data/elections/PR2012_T2.csv", sep="\t", colClasses = c(rep("character", 2), rep("numeric", 16)), head = TRUE, stringsAsFactors = TRUE, dec =",")

variables_elect <- c("PR2012_T1_Abst_prct_insc", "PR2012_T1_prct_insc_LE.PEN", "PR2012_T1_prct_insc_SARKOZY", "PR2012_T1_prct_insc_MÉLENCHON", "PR2012_T1_prct_insc_BAYROU", "PR2012_T1_prct_insc_HOLLANDE")

colPR2012_T1 <- colnames(PR2012_T1)
# colPR2012_T2 <- colnames(PR2012_T2)

# On liste les variables présentes seulement dans budgets (celles qu'il va falloir y prélever)
colgarder2011_PR2012_T1 <- setdiff(colPR2012_T1, coldf2011)
# colgarder2011_PR2012_T2 <- setdiff(colPR2012_T2, coldf2011)

# On fusionne les données de data_prep et des budgets, en ne gardant pour ces dernières que les données précédemment identifiées
df2011 <- merge(df2011,
                     PR2012_T1[, c("CODGEO", colgarder2011_PR2012_T1)],
                     by = "CODGEO", all.x = TRUE, all.y = TRUE)
# Suppression des données inutiles (déjà importées ou seulement de transition)
# df2011 <- merge(df2011,
#                      PR2012_T2[, c("CODGEO", colgarder2011_PR2012_T2)],
#                      by = "CODGEO", all.x = TRUE, all.y = TRUE)
rm(PR2012_T1, colPR2012_T1, coldf2011, colgarder2011_PR2012_T1)

3.5 Données sur les RPI

Données pas utilisables en l’état (sont à la géométrie 2024).

load("data/RPI_par_communes_2024.Rdata")
df2011 <- merge(df2011, RPI_par_communes,
                     by = "CODGEO", all.x = TRUE)
table(df2011$RPI)
summary(is.na(df2011$RPI))
df2011$RPI[is.na(df2011$RPI)] <- "Sans_ecoles" # On indique comme étant "Sans_ecoles" les communes pour lesquelles aucune école n'était renseignée dans l'annuaire de l'Éducation nationale.
table(df2011$RPI)
summary(is.na(df2011$RPI))

3.6 Données sur la BPE (Banque Permanente des Équipements)

load("data/BPE_2014.Rdata")
# On remplace l'absence de valeurs par des zéros
# BPE_communes[is.na(BPE_communes)] <- 0
rm(BPE_communes)
BPE_communes_decode[is.na(BPE_communes_decode)] <- 0

# Pour alléger les donées, on se concentre sur quelques équipements uniquements
BPE_communes_decode <- BPE_communes_decode[, c("CODGEO", "Coiffure", "Supermarché", "École élémentaire", "Médecin omnipraticien", "Taxi", "Salles multisports (gymnase)", "Infirmier", "Pharmacie")]

# Il faut ensuite recoder les variables car on veut faire des catégories
equipement <- "Médecin omnipraticien"
for (equipement in colnames(BPE_communes_decode)[2:ncol(BPE_communes_decode)]) {
  seuils <- quantile(BPE_communes_decode[, equipement], probs = seq(0, 1, 0.05), na.rm = TRUE) # Définition des seuils par déciles
  seuils <- unique(seuils) # On ne prend que des valeurs uniques
  seuils[length(seuils)] <- seuils[length(seuils)]+1 # On rajoute 1 à la dernière valeur pour que les catégories définies ensuite englobent toutes les communes
  BPE_communes_decode[, paste0(equipement, "_quantile")] <- cut(BPE_communes_decode[, equipement], breaks = seuils, right = FALSE) # right : pour ouverture de la coupure vers la gauche et non vers la droite comme par défait
    round(100*prop.table(table(BPE_communes_decode[, equipement]),margin=),2)
    round(100*prop.table(table(BPE_communes_decode[, paste0(equipement, "_quantile")]),margin=),2)

  levels(BPE_communes_decode[, paste0(equipement, "_quantile")])
  
}

df2011 <- merge(df2011, BPE_communes_decode,
                     by = "CODGEO", all.x = TRUE)


# table(BPE_communes_decode$Supermarché)
# round(100*prop.table(table(BPE_communes_decode$Supermarché),margin=),2)

4 Préparation des données

4.1 Sélection études de cas

Section qui sert, au départ, à tester le code sur de petits ensembles.

À l’avenir, pourra servir à définir des boucles

EdCs <- c("EdC_France", "EdC_Savoies", "EdC_Ouest", "EdC_Normandie", "EdC_49", "EdC_RALP_partiel")

df2011$EdC_France <- "OUI" # France entière

df2011$EdC_Savoies[df2011$CODE_DEPT %in% c("73", "74")] <- "OUI" # Départements de Savoie et Haute-Savoie
df2011[, c("EdC_Savoies")] [is.na(df2011[, c("EdC_Savoies")])] <- "NON"

df2011$EdC_Ouest[df2011$REG %in% c("23", "25", "53", "52")] <- "OUI" # Haute-Normandie, Basse Normandie, Bretagne, Pays-de-la-Loire
df2011[, c("EdC_Ouest")] [is.na(df2011[, c("EdC_Ouest")])] <- "NON"

df2011$EdC_Normandie[df2011$REG %in% c("23", "25")] <- "OUI" # Haute et Basse Normandie
df2011[, c("EdC_Normandie")] [is.na(df2011[, c("EdC_Normandie")])] <- "NON"

df2011$EdC_49[df2011$CODE_DEPT %in% c("49")] <- "OUI" # Département du Maine-et-Loire
df2011[, c("EdC_49")] [is.na(df2011[, c("EdC_49")])] <- "NON"

df2011$EdC_RALP_partiel[df2011$CODE_DEPT %in% c("69", "73", "74", "38", "01")] <- "OUI" # Départements de Savoie, Haute-Savoie, Isère, Rhône et Ain
df2011[, c("EdC_RALP_partiel")] [is.na(df2011[, c("EdC_RALP_partiel")])] <- "NON"

EdC <- EdCs[1] # France
# EdC <- EdCs[3] # Ouest
# EdC <- EdCs[4] # Normandie
# EdC <- EdCs[5] # Maine-et-Loire

df_reg <- subset(df2011, df2011[EdC] == "OUI")
df_reg_Cfus <- subset(df_reg, df_reg$COM_NOUV == "OUI")

# On peut choisir de se limiter aux communes françaises ayant une population ne dépassant pas celle de la plus peuplée des communes nouvelles
df_reg <- subset(df_reg, df_reg$P09_POP <= max(df_reg_Cfus$P09_POP))

# Si on veut ne sélectionner que certains départements 
# departementsEdC <- as.data.frame(table(df_reg$CODE_DEPT))[1]
# dep <- merge(dep, departementsEdC, by.x = "CODE_DEPT", by.y = "Var1", all.x = FALSE)

# On ajoute les géométries
geomfus2011 <- merge(geomfus2011, df_reg[, c("CODGEO", "LIBGEO")], by = "CODGEO", all.x = FALSE, all.y = TRUE)
df_reg <- merge(geom2011, df_reg, by = "CODGEO", all.y = TRUE)



plot(df_reg$geom, main = paste("Espace étudié :", EdC))

4.2 Adaptation variables pourcentages [optionnel]

NB : Il y a un problème si on utilise les variables budgétaires en valeurs absolues ou les variables en pourcentages notées 100% = 100. Cela fonctionne plus correctement si on a les pourcentages notées en 100% = 1.

Deux possibilités : diviser les pourcentages manuellement ou utiliser le fichier ‘refdata_budgets_prct_en_ratio.Rdata’. On choisit la première solution.

Pas vraiment nécessaire si passage des variables quantitatives en variables qualitatives (quantiles, cf. section suivante).

# Pour multiplier/diviser
# Toutes les variables en pourcentages qui sont des variables budgétaires
variables_budg_prct <- colnames(df_reg[which(stringr::str_detect(colnames(df_reg), "_prct") == TRUE)])
variables_budg_prct <- variables_budg_prct[variables_budg_prct!= "geometry"] # On supprime la variable "geometry" qui est intégrée sans que cela soit voulu
df_reg [variables_budg_prct] <- df_reg[variables_budg_prct]/100
# test <- df_reg[variables_budg_prct]

variables_RT <- colnames(df_reg[which(stringr::str_detect(colnames(df_reg), "_RT") == TRUE)])
variables_RT <- variables_RT[variables_RT!= "geometry"] # On supprime la variable "geometry" qui est intégrée sans que cela soit voulu
df_reg [variables_RT] <- df_reg[variables_RT]/100
# test <- df_reg[variables_RT]

# Sans que je ne comprenne bien pourquoi, les géométries sont également touchées, si variable non retirée avant, on rétablit donc ces dernières
# df_reg$geometry <- df_reg$geometry * 10000

4.3 Transformations des variables quantitatives en variables qualitatives

Le principe d’une régression logistique est de définir la probabilité d’un évènement (dans mon cas, fusionner ou non) en fonction du changement d’une variable. Dans le cas d’une variable quantitative, la question est quelle est le pourcentage de chance que l’évènement survienne davantage lorsqu’on augmente la variable quantitative de 1.

Dans le cas d’une variable qualitative, il s’agit du pourcentage de chance que l’évènement survienne davantage lorsqu’on change de catégorie. Il faut donc faire attention à mettre la variable en facteur, avec la première valeur qui doit être celle de référence (c’est à partir de celle-là que les comparaisons vont être faites, c’est donc important car c’est cela qui détermine le discours).

# variable <- "C09_ACTOCC_OUT_RT"
variables_quanti_quali <- c(variables_budg_prct, variables_RT, variables_elect, "superficie")
# On retire certaines variables qui ne peuvent être divisées par tertiles car trop de valeurs nulles.
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DSU_tot_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DNP_princip_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DNP_sortie_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DNP_major_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DNP_tot_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DSR_bourg_centre_global_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DSR_cible_global_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="C09_EMPLT_INDUS_RT"]

df_sans_geom <- df_reg # Création d'un jeu de données pour le travail
st_geometry(df_sans_geom) <- NULL # Suppression des géométries
# summary(df_sans_geom[, variable])
variable <- "superficie"
df_sans_geom$superficie <- as.numeric(df_sans_geom$superficie)
# class(df_sans_geom$superficie)

for (variable in variables_quanti_quali) {

Y <- df_sans_geom[, variable] # On se focalise sur une variable en particulier

# Calcul moyenne groupes
moy <- as.data.frame(tapply(df_sans_geom[, variable], df_sans_geom$COM_NOUV, mean, na.rm = TRUE))
moy$COM_NOUV <- row.names(moy)
colnames(moy) <- c("moyenne", "COM_NOUV") 

# Calcul moyenne groupes
med <- as.data.frame(tapply(df_sans_geom[, variable], df_sans_geom$COM_NOUV, median, na.rm = TRUE))
med$COM_NOUV <- row.names(moy)
colnames(med) <- c("median", "COM_NOUV") 

# On veut observer graphiquement la distribution statistique
histo <- ggplot(df_sans_geom)+
  geom_histogram(aes(x=df_sans_geom[, variable], color = COM_NOUV), bins=50)+
  labs(title = paste("Variable étudiée :\n", variable), subtitle = "Pointillés : moyennes.\nLigne pleine : médianes.")+
  # geom_vline(aes(xintercept=median(df_sans_geom[, variable], na.rm = T)),      color="black", linetype="dashed" )+
  geom_vline(data=moy, aes(xintercept=moyenne, color=COM_NOUV),             linetype="dashed") +
  geom_vline(data=med, aes(xintercept=median, color=COM_NOUV)) +
  ylab("Count")+theme_light()
print(histo)

# On découpe en fonction du quantile souhaité puis on modifie l'ordre de niveaux de facteur pour que la variable centrale soit la première

coupures <- quantile(Y, probs=seq(0, 1, 0.2), na.rm = TRUE)

if (any(duplicated(coupures)) == FALSE) { # Si c'est possible, on utilise les quintiles
  Y <- cut(Y,breaks=c(quantile(Y, probs=seq(0, 1, 0.2), na.rm = TRUE))) # Pour découpage en quintiles
levels(Y)<-c("Q1","Q2", "Q3", "Q4", "Q5")
Y <- relevel (Y, "Q3","Q1", "Q2", "Q4", "Q5")
}else{# Sinon, on utilise les tertiles
Y <- cut(Y,breaks=c(quantile(Y, probs=seq(0, 1, 1/3), na.rm = TRUE))) # Pour découpage en tertiles
levels(Y)<-c("Tertile_inf","Tertile_med", "Tertile_sup")
Y <- relevel (Y, "Tertile_med", "Tertile_inf", "Tertile_sup")
}


# Y<-cut(Y,breaks=c(quantile(Y))) # Pour découpage en quartiles
# levels(Y)<-c("Q1","Q2", "Q3", "Q4")




summary(Y)


X <- df_sans_geom$COM_NOUV
summary(X)

tabcont<-table(X,Y)
tabcont # En valeur absolue
# round(100*prop.table(tabcont,margin=1),1) # Pourcentages, le total se fait par lignes
# round(100*prop.table(tabcont,margin=),1) # Pourcentages, le total se fait sur l'ensemble de la population
# round(100*prop.table(tabcont,margin=2),1) # Pourcentages, le total se fait par colonnes

# On verse les nouvelles données au data frame
df_reg[, paste0(variable, "_quali")] <- Y

}

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 52 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 52 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 52 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 101 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 126 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 958 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

rm(tabcont, X, Y, df_sans_geom, variable, histo)

4.4 Modifications de quelques autres variables

Là encore, pour les transformer en variables qualitatives.

# On réordonne le facteur CATAEU2010
df_reg$CATAEU2010 <- as.factor(df_reg$CATAEU2010)
freq(df_reg$CATAEU2010, sort = "dec")

##         n    % val%
## 112 12180 33.7 33.7
## 400  7233 20.0 20.0
## 300  6975 19.3 19.3
## 120  3968 11.0 11.0
## 111  3100  8.6  8.6
## 221   854  2.4  2.4
## 212   801  2.2  2.2
## 222   558  1.5  1.5
## 211   439  1.2  1.2

row.names(freq(df_reg$CATAEU2010, sort = "dec"))

## [1] "112" "400" "300" "120" "111" "221" "212" "222" "211"

# Vérifier que la modalité qui se trouve en premier est bien la plus fréquente.
# Sinon, possibilité d'utiliser le code suivant :
# df_reg$CATAEU2010 <- relevel (df_reg$CATAEU2010, "112", "111", "120", "211", "212", "221", "222", "300", "400")
# df_reg$CATAEU2010 <- relevel (df_reg$CATAEU2010, "112", "400", "300", "120", "111", "221", "212", "222", "211")
ordre <- row.names(freq(df_reg$CATAEU2010, sort = "dec"))
df_reg$CATAEU2010 <- factor (df_reg$CATAEU2010, ordre)


df_reg$CODE_DEPT <- as.factor(df_reg$CODE_DEPT)
freq(df_reg$CODE_DEPT, sort = "dec")

##      n   % val%
## 62 894 2.5  2.5
## 02 815 2.3  2.3
## 80 781 2.2  2.2
## 76 743 2.1  2.1
## 57 729 2.0  2.0
## 14 705 2.0  2.0
## 21 705 2.0  2.0
## 60 692 1.9  1.9
## 27 675 1.9  1.9
## 59 645 1.8  1.8
## 51 619 1.7  1.7
## 50 601 1.7  1.7
## 25 593 1.6  1.6
## 54 593 1.6  1.6
## 31 588 1.6  1.6
## 71 573 1.6  1.6
## 24 557 1.5  1.5
## 64 546 1.5  1.5
## 70 545 1.5  1.5
## 39 544 1.5  1.5
## 33 539 1.5  1.5
## 38 532 1.5  1.5
## 67 526 1.5  1.5
## 88 515 1.4  1.4
## 77 513 1.4  1.4
## 61 505 1.4  1.4
## 55 500 1.4  1.4
## 65 474 1.3  1.3
## 17 471 1.3  1.3
## 63 469 1.3  1.3
## 08 463 1.3  1.3
## 32 463 1.3  1.3
## 89 455 1.3  1.3
##  [ reached 'max' / getOption("max.print") -- omitted 60 rows ]

ordre <- row.names(freq(df_reg$CODE_DEPT, sort = "dec"))
df_reg$CODE_DEPT <- factor (df_reg$CODE_DEPT, ordre)

df_reg$REG <- as.factor(df_reg$REG)
freq(df_reg$REG, sort = "dec")

##       n   % val%
## 73 3018 8.4  8.4
## 82 2872 8.0  8.0
## 41 2337 6.5  6.5
## 72 2292 6.3  6.3
## 22 2288 6.3  6.3
## 26 2045 5.7  5.7
## 21 1947 5.4  5.4
## 24 1839 5.1  5.1
## 25 1811 5.0  5.0
## 43 1784 4.9  4.9
## 91 1541 4.3  4.3
## 31 1539 4.3  4.3
## 52 1497 4.1  4.1
## 54 1459 4.0  4.0
## 23 1418 3.9  3.9
## 83 1309 3.6  3.6
## 53 1265 3.5  3.5
## 11 1247 3.5  3.5
## 93  953 2.6  2.6
## 42  901 2.5  2.5
## 74  746 2.1  2.1

row.names(freq(df_reg$REG, sort = "dec"))

##  [1] "73" "82" "41" "72" "22" "26" "21" "24" "25" "43" "91" "31" "52" "54" "23"
## [16] "83" "53" "11" "93" "42" "74"

ordre <- row.names(freq(df_reg$REG, sort = "dec"))
df_reg$REG <- factor (df_reg$REG, ordre)


df_reg$P09_POP <- as.numeric(df_reg$P09_POP)
df_reg$P09_POP_quali <- cut(df_reg$P09_POP, breaks=c(quantile(df_reg$P09_POP, probs=seq(0, 1, 1/5), na.rm = TRUE)))
levels(df_reg$P09_POP_quali)<-c("Q1","Q2", "Q3", "Q4", "Q5")
df_reg$P09_POP_quali <- relevel (df_reg$P09_POP_quali, "Q3","Q1", "Q2", "Q4", "Q5")

# RPI pas utilisés pour l'instant
# df_reg$RPI <- as.factor(df_reg$RPI)
# freq(df_reg$RPI, sort = "dec")
# row.names(freq(df_reg$RPI, sort = "dec"))
# ordre <- row.names(freq(df_reg$RPI, sort = "dec"))
# ordre <- c("Ecoles_sans_RPI", "Sans_ecoles", "RPI", "Avec_et_sans_RPI")
# df_reg$RPI <- factor (df_reg$RPI, ordre)

5 Étude auto-corrélation et colinéarité

Pour éviter la colinéarité entre les variables explicatives, deux méthodes [feuillet2019] :

faire une matrice de corrélation et supprimer les variables dont les corrélations sont supérieures à 0,7 en valeur absolue [tenenhaus2007] ;
calculer la Variance Inflation Factor et supprimer les variables avec une VIF supérieure à 3 (voire supérieur à 2, même référence biblio).

5.1 Réalisation d’une matrice de corrélation

Cf. https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html

# Pour sélectionner certaines variables spécifiquement
# variables_correlations <- c("fctva_prct", "dgf_prct", "DSU_tot_prct", "DNP_tot_prct", "DSR_pereq_global_prct")
# variables_correlations <- c("fctva", "dgf", "DSU_tot", "DNP_tot", "DSR_pereq_global")
# variables_correlations <- vartests

# Toutes les variables en pourcentages NB : ajout du symbole $ pour indiquer que c'est la fin de la chaîne de caractère
variables_correlations <- c(colnames(df_reg[which(stringr::str_detect(colnames(df_reg), "_prct$") == TRUE)]),
                            colnames(df_reg[which(stringr::str_detect(colnames(df_reg), "_RT$") == TRUE)]))
variables_correlations <- c(colnames(df_reg[which(stringr::str_detect(colnames(df_reg), "_prct$") == TRUE)]),
                            colnames(df_reg[which(stringr::str_detect(colnames(df_reg), "_RT$") == TRUE)]), variables_elect)

variables_correlations <- variables_correlations[variables_correlations!= "geometry"] # On supprime la variable "geometry" si elle est intégrée sans que cela soit voulu

# variables_correlations <- c("P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "P09_POP3044Y_RT", "C09_ACTOCC_OUT_RT", "P09_CHOM1564_RT", "dgf_prct", "perso_prct", "charge_prct")

# Définition du seuil d'auto-corrélation
seuil <- 0.7

# Juste pour les communes fusionnées
tmp <- subset(df2011, df2011$COM_NOUV == "OUI")
tmp <- na.omit(tmp [, variables_correlations])

tmp2 <- cor(tmp, method = "pearson")
res1 <- cor.mtest(tmp, conf.level = .99)

col <- c("#a50026","#d73027","#f46d43","#fdae61","#fee090","#e0f3f8","#abd9e9","#74add1","#4575b4","#313695")

corrplot::corrplot(tmp2, p.mat = res1$p, sig.level = .2, order = "hclust", tl.cex = 0.8, tl.col = "black") # Obligé d'appeler le package directement car sinon conflit avec une autre fonction du package `arm`.

corrplot::corrplot(tmp2, method="color", order = "hclust", # cl.lim = c(-1,1), # cl.lim pose problème donc supprimé
         col = col,
         tl.col = rgb(0,0,0),
         number.cex = 0.4,
         addCoef.col = rgb(0,0,0, alpha =0.6),
         addgrid.col = "white",
         title = "Analyse des corrélations",
         tl.cex=0.8,
         cl.cex = 0.6,
         mar=c(0,0,4,0))

tmp2[which(tmp2 == 1)] <- 0 # On remplace par 0 l'autocorrélations des variables avec elles-mêmes pour faciliter la lecture ensuite
# which(abs(tmp2) > 0.7, arr.ind=TRUE) # On identifie les valeurs dont la valeur absolue serait supérieure à 0,7

for(i in 1:nrow(tmp2)) {
   for(j in 1:ncol(tmp2))
   {if(tmp2[i,j]>seuil)
   {print(paste("Variables auto-corrélées à plus de", seuil, ":", rownames(tmp2)[i], "et", colnames(tmp2)[j], ":", tmp2[i,j]))
     }}}

## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et charge_prct : 0.732689123376272"
## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et res1_prct : 0.734728134238698"
## [1] "Variables auto-corrélées à plus de 0.7 : impo1_prct et bfb_prct : 0.934797223269612"
## [1] "Variables auto-corrélées à plus de 0.7 : impo1_prct et pfb_prct : 0.944226808552536"
## [1] "Variables auto-corrélées à plus de 0.7 : impo1_prct et PotFin4taxes_prct : 0.723731388000992"
## [1] "Variables auto-corrélées à plus de 0.7 : impo1_prct et P11_POT_FIN_RT : 0.823753548373962"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DotationNmoins1PerimN_prct : 0.958902817093769"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et Dotation_forfaitaireN_prct : 0.96347828612023"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DSR_pereq_global_prct : 0.826044221470363"
## [1] "Variables auto-corrélées à plus de 0.7 : charge_prct et prod_prct : 0.732689123376272"
## [1] "Variables auto-corrélées à plus de 0.7 : res1_prct et prod_prct : 0.734728134238698"
## [1] "Variables auto-corrélées à plus de 0.7 : recinv_prct et res2_prct : 0.741258939982701"
## [1] "Variables auto-corrélées à plus de 0.7 : depinv_prct et equip_prct : 0.887426789236825"
## [1] "Variables auto-corrélées à plus de 0.7 : equip_prct et depinv_prct : 0.887426789236825"
## [1] "Variables auto-corrélées à plus de 0.7 : remb_prct et annu_prct : 0.97130987424642"
## [1] "Variables auto-corrélées à plus de 0.7 : res2_prct et recinv_prct : 0.741258939982701"
## [1] "Variables auto-corrélées à plus de 0.7 : annu_prct et remb_prct : 0.97130987424642"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et impo1_prct : 0.934797223269612"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et pfb_prct : 0.987194884577028"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et PotFin4taxes_prct : 0.722401553651033"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et P11_POT_FIN_RT : 0.87843217078898"
## [1] "Variables auto-corrélées à plus de 0.7 : bfnb_prct et pfnb_prct : 0.736196383226621"
## [1] "Variables auto-corrélées à plus de 0.7 : pfb_prct et impo1_prct : 0.944226808552536"
## [1] "Variables auto-corrélées à plus de 0.7 : pfb_prct et bfb_prct : 0.987194884577028"
## [1] "Variables auto-corrélées à plus de 0.7 : pfb_prct et P11_POT_FIN_RT : 0.884623222219747"
## [1] "Variables auto-corrélées à plus de 0.7 : pfnb_prct et bfnb_prct : 0.736196383226621"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et dgf_prct : 0.958902817093769"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et Dotation_forfaitaireN_prct : 0.996833840765774"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et DSR_pereq_global_prct : 0.804191406661557"
## [1] "Variables auto-corrélées à plus de 0.7 : PotFin4taxes_prct et impo1_prct : 0.723731388000992"
## [1] "Variables auto-corrélées à plus de 0.7 : PotFin4taxes_prct et bfb_prct : 0.722401553651033"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et dgf_prct : 0.96347828612023"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DotationNmoins1PerimN_prct : 0.996833840765774"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DSR_pereq_global_prct : 0.812335469877741"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_major_prct : 0.780156401039296"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_tot_prct : 0.97878473516058"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_princip_prct : 0.780156401039296"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_tot_prct : 0.890785442109139"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_princip_prct : 0.97878473516058"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_major_prct : 0.890785442109139"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et dgf_prct : 0.826044221470363"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et DotationNmoins1PerimN_prct : 0.804191406661557"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et Dotation_forfaitaireN_prct : 0.812335469877741"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_RETR1564_RT et P09_POP6074Y_RT : 0.742995991320959"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_ACT1564_Agr_RT et C09_EMPLT_AGRI_RT : 0.722572663768533"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_EMPLT_AGRI_RT et C09_ACT1564_Agr_RT : 0.722572663768533"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_POP6074Y_RT et P09_RETR1564_RT : 0.742995991320959"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_POT_FIN_RT et impo1_prct : 0.823753548373962"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_POT_FIN_RT et bfb_prct : 0.87843217078898"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_POT_FIN_RT et pfb_prct : 0.884623222219747"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_IMP_NET_RT : 0.861018323277089"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_FoyFisc_Imp_RT : 0.738703335772406"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_IMP_NET_RT et P11_Rev_Fisc_RT : 0.861018323277089"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_FoyFisc_Imp_RT et P11_Rev_Fisc_RT : 0.738703335772406"

# Communes n'ayant pas fusionné
tmp <- subset(df2011, df2011$COM_NOUV == "NON")
tmp <- na.omit(tmp [, variables_correlations])

tmp2 <- cor(tmp, method = "pearson")
res1 <- cor.mtest(tmp, method = "pearson")

corrplot::corrplot(tmp2, p.mat = res1$p, sig.level = .2, order = "hclust",
         tl.cex = 0.8, tl.col = "black")

col <- c("#a50026","#d73027","#f46d43","#fdae61","#fee090","#e0f3f8","#abd9e9","#74add1","#4575b4","#313695")

corrplot::corrplot(tmp2, method="color", order = "hclust", # cl.lim = c(-1,1), # cl.lim pose problème donc supprimé 
         col = col,
         tl.col = rgb(0,0,0),
         number.cex = 0.4,
         addCoef.col = rgb(0,0,0, alpha =0.6),
         addgrid.col = "white",
         title = "Analyse des corrélations",
         tl.cex=0.8,
         cl.cex = 0.6,
         mar=c(0,0,4,0))

tmp2[which(tmp2 == 1)] <- 0 # On remplace par 0 l'autocorrélations des variables avec elles-mêmes pour faciliter la lecture ensuite
# which(abs(tmp2) > 0.7, arr.ind=TRUE) # On identifie les valeurs dont la valeur absolue serait supérieure à 0,7

for(i in 1:nrow(tmp2)) {
   for(j in 1:ncol(tmp2))
   {if(tmp2[i,j]>seuil)
   {print(paste("Variables auto-corrélées à plus de", seuil, ":", rownames(tmp2)[i], "et", colnames(tmp2)[j], ":", tmp2[i,j]))
     }}}

## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et charge_prct : 0.724234792513559"
## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et res1_prct : 0.732830571012828"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DotationNmoins1PerimN_prct : 0.945250511144103"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et Dotation_forfaitaireN_prct : 0.953487917133836"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DSR_pereq_global_prct : 0.79232599936631"
## [1] "Variables auto-corrélées à plus de 0.7 : charge_prct et prod_prct : 0.724234792513559"
## [1] "Variables auto-corrélées à plus de 0.7 : res1_prct et prod_prct : 0.732830571012828"
## [1] "Variables auto-corrélées à plus de 0.7 : recinv_prct et res2_prct : 0.723036360221147"
## [1] "Variables auto-corrélées à plus de 0.7 : depinv_prct et equip_prct : 0.902821742043709"
## [1] "Variables auto-corrélées à plus de 0.7 : equip_prct et depinv_prct : 0.902821742043709"
## [1] "Variables auto-corrélées à plus de 0.7 : remb_prct et annu_prct : 0.971637705415465"
## [1] "Variables auto-corrélées à plus de 0.7 : res2_prct et recinv_prct : 0.723036360221147"
## [1] "Variables auto-corrélées à plus de 0.7 : annu_prct et remb_prct : 0.971637705415465"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et PotFin4taxes_prct : 0.782069994452378"
## [1] "Variables auto-corrélées à plus de 0.7 : bfnb_prct et pfnb_prct : 0.733021436204166"
## [1] "Variables auto-corrélées à plus de 0.7 : pfnb_prct et bfnb_prct : 0.733021436204166"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et dgf_prct : 0.945250511144103"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et Dotation_forfaitaireN_prct : 0.991260468887141"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et DSR_pereq_global_prct : 0.771690289546452"
## [1] "Variables auto-corrélées à plus de 0.7 : PotFin4taxes_prct et bfb_prct : 0.782069994452378"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et dgf_prct : 0.953487917133836"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DotationNmoins1PerimN_prct : 0.991260468887141"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DSR_pereq_global_prct : 0.777669035129182"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_major_prct : 0.79425679995476"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_tot_prct : 0.979844133560585"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_princip_prct : 0.79425679995476"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_tot_prct : 0.897902584225198"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_princip_prct : 0.979844133560585"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_major_prct : 0.897902584225198"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et dgf_prct : 0.79232599936631"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et DotationNmoins1PerimN_prct : 0.771690289546452"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et Dotation_forfaitaireN_prct : 0.777669035129182"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_RETR1564_RT et P09_POP6074Y_RT : 0.740842384941785"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_ACT1564_Agr_RT et C09_EMPLT_AGRI_RT : 0.713423404098747"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_EMPLT_AGRI_RT et C09_ACT1564_Agr_RT : 0.713423404098747"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_POP6074Y_RT et P09_RETR1564_RT : 0.740842384941785"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_IMP_NET_RT : 0.880125493065463"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_FoyFisc_Imp_RT : 0.734113201668045"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_IMP_NET_RT et P11_Rev_Fisc_RT : 0.880125493065463"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_FoyFisc_Imp_RT et P11_Rev_Fisc_RT : 0.734113201668045"

# Ensemble des communes françaises
tmp <- df2011
tmp <- na.omit(tmp [, variables_correlations])

tmp2 <- cor(tmp, method = "pearson")
res1 <- cor.mtest(tmp, method = "pearson")

corrplot::corrplot(tmp2, p.mat = res1$p, sig.level = .2, order = "hclust",
         tl.cex = 0.8, tl.col = "black")

col <- c("#a50026","#d73027","#f46d43","#fdae61","#fee090","#e0f3f8","#abd9e9","#74add1","#4575b4","#313695")

corrplot::corrplot(tmp2, method="color", order = "hclust", # cl.lim = c(-1,1), # cl.lim pose problème donc supprimé 
         col = col,
         tl.col = rgb(0,0,0),
         number.cex = 0.4,
         addCoef.col = rgb(0,0,0, alpha =0.6),
         addgrid.col = "white",
         title = "Analyse des corrélations",
         tl.cex=0.8,
         cl.cex = 0.6,
         mar=c(0,0,4,0))

tmp2[which(tmp2 == 1)] <- 0 # On remplace par 0 l'autocorrélations des variables avec elles-mêmes pour faciliter la lecture ensuite
# which(abs(tmp2) > 0.7, arr.ind=TRUE) # On identifie les valeurs dont la valeur absolue serait supérieure à 0,7

for(i in 1:nrow(tmp2)) {
   for(j in 1:ncol(tmp2))
   {if(tmp2[i,j]>seuil)
   {print(paste("Variables auto-corrélées à plus de", seuil, ":", rownames(tmp2)[i], "et", colnames(tmp2)[j], ":", tmp2[i,j]))
     }}}

## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et charge_prct : 0.724818656916716"
## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et res1_prct : 0.732938988586129"
## [1] "Variables auto-corrélées à plus de 0.7 : impo1_prct et pfb_prct : 0.719361456090989"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DotationNmoins1PerimN_prct : 0.946561903962206"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et Dotation_forfaitaireN_prct : 0.954457043664577"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DSR_pereq_global_prct : 0.795902189919594"
## [1] "Variables auto-corrélées à plus de 0.7 : charge_prct et prod_prct : 0.724818656916716"
## [1] "Variables auto-corrélées à plus de 0.7 : res1_prct et prod_prct : 0.732938988586129"
## [1] "Variables auto-corrélées à plus de 0.7 : recinv_prct et res2_prct : 0.724291657635224"
## [1] "Variables auto-corrélées à plus de 0.7 : depinv_prct et equip_prct : 0.901768326955515"
## [1] "Variables auto-corrélées à plus de 0.7 : equip_prct et depinv_prct : 0.901768326955515"
## [1] "Variables auto-corrélées à plus de 0.7 : remb_prct et annu_prct : 0.971615494510789"
## [1] "Variables auto-corrélées à plus de 0.7 : res2_prct et recinv_prct : 0.724291657635224"
## [1] "Variables auto-corrélées à plus de 0.7 : annu_prct et remb_prct : 0.971615494510789"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et pfb_prct : 0.878040662678628"
## [1] "Variables auto-corrélées à plus de 0.7 : bfnb_prct et pfnb_prct : 0.734066280198133"
## [1] "Variables auto-corrélées à plus de 0.7 : pfb_prct et impo1_prct : 0.719361456090989"
## [1] "Variables auto-corrélées à plus de 0.7 : pfb_prct et bfb_prct : 0.878040662678628"
## [1] "Variables auto-corrélées à plus de 0.7 : pfnb_prct et bfnb_prct : 0.734066280198133"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et dgf_prct : 0.946561903962206"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et Dotation_forfaitaireN_prct : 0.991748574148787"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et DSR_pereq_global_prct : 0.775196322883291"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et dgf_prct : 0.954457043664577"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DotationNmoins1PerimN_prct : 0.991748574148787"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DSR_pereq_global_prct : 0.781293401890802"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_major_prct : 0.793043813747828"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_tot_prct : 0.979764635015164"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_princip_prct : 0.793043813747828"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_tot_prct : 0.897271624224915"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_princip_prct : 0.979764635015164"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_major_prct : 0.897271624224915"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et dgf_prct : 0.795902189919594"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et DotationNmoins1PerimN_prct : 0.775196322883291"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et Dotation_forfaitaireN_prct : 0.781293401890802"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_RETR1564_RT et P09_POP6074Y_RT : 0.740932400425599"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_ACT1564_Agr_RT et C09_EMPLT_AGRI_RT : 0.714166544320671"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_EMPLT_AGRI_RT et C09_ACT1564_Agr_RT : 0.714166544320671"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_POP6074Y_RT et P09_RETR1564_RT : 0.740932400425599"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_IMP_NET_RT : 0.879492075658382"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_FoyFisc_Imp_RT : 0.734321656244132"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_IMP_NET_RT et P11_Rev_Fisc_RT : 0.879492075658382"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_FoyFisc_Imp_RT et P11_Rev_Fisc_RT : 0.734321656244132"

Les variables envisagées ne sont pas trop auto-corrélées entre elles si on se réfère aux seuils préconisés dans la littérature.

Si jamais des variables sont trop autocorrélées, on peut les supprimer (c’est ce qui est fait dans la section suivante) puis relancer une étude de l’auto-corrélation.

variables_correlations <- variables_correlations[variables_correlations!="impo1_prct"]
variables_correlations <- variables_correlations[variables_correlations!="Pop_INSEE_prct"]
variables_correlations <- variables_correlations[variables_correlations!="Pop_DGF_prct"]
variables_correlations <- variables_correlations[variables_correlations!="Dotation_forfaitaireN_prct"]
variables_correlations <- variables_correlations[variables_correlations!="pfb_prct"]
variables_correlations <- variables_correlations[variables_correlations!="DNP_princip_prct"]
variables_correlations <- variables_correlations[variables_correlations!="DNP_major_prct"]
variables_correlations <- variables_correlations[variables_correlations!="DNP_tot_prct"]
variables_correlations <- variables_correlations[variables_correlations!="DSR_pereq_global_prct"]
variables_correlations <- variables_correlations[variables_correlations!="P11_POT_FIN_RT"]
variables_correlations <- variables_correlations[variables_correlations!="P11_Rev_Fisc_RT"]
variables_correlations <- variables_correlations[variables_correlations!="P11_IMP_NET_RT"]
variables_correlations <- variables_correlations[variables_correlations!="res1_prct"]
variables_correlations <- variables_correlations[variables_correlations!="depinv_prct"]
variables_correlations <- variables_correlations[variables_correlations!="tth_prct"]
variables_correlations <- variables_correlations[variables_correlations!="prod_prct"]
variables_correlations <- variables_correlations[variables_correlations!="res2_prct"]
variables_correlations <- variables_correlations[variables_correlations!="remb_prct"]
variables_correlations <- variables_correlations[variables_correlations!="pfnb_prct"]

rm(tmp, tmp2, res1, col, i, j, seuil)

## Warning in rm(tmp, tmp2, res1, col, i, j, seuil): objet 'tmp' introuvable

## Warning in rm(tmp, tmp2, res1, col, i, j, seuil): objet 'tmp2' introuvable

## Warning in rm(tmp, tmp2, res1, col, i, j, seuil): objet 'res1' introuvable

## Warning in rm(tmp, tmp2, res1, col, i, j, seuil): objet 'col' introuvable

## Warning in rm(tmp, tmp2, res1, col, i, j, seuil): objet 'i' introuvable

## Warning in rm(tmp, tmp2, res1, col, i, j, seuil): objet 'j' introuvable

5.2 Calcul de la Variance Inflaction Factor (VIF)

D’après la littérature [tenenhaus2007] , il faut calculer la Variance Inflation Factor (VIF) et supprimer les variables avec une VIF supérieure à 3 (voire supérieur à 2).

Pour la mise en œuvre, cf. cette page : https://larmarange.github.io/guide-R/analyses/multicolinearite.html

Ces éléments seront donc présentés lors de la mise en œuvre de la régression.

6 Comparaison des valeurs des communes fusionnantes vis-à-vis des valeurs des communes inchangées

# Calculs des valeurs moyennes (pour comparaison)
# On identifie les variables en stock
variables_stock <- stringr::str_replace(variables_RT, "Y_RT", "")
variables_stock <- stringr::str_replace(variables_stock, "_RT", "")
variables_stock <- variables_stock[-which(variables_stock == "C09_EMP_CONC")]
variables_stock <- variables_stock[1:(length(variables_stock)-5)]

# Import données totales
# load("data/refdata.Rdata")
# rm(df_new)
# Sélection données utiles
pourmoyennesCAH <- na.omit(subset(df2011, COM_NOUV == "OUI"))
moyennesCAH <- data.frame()
i <- variables_stock[2]
ratio <- as.data.frame(read_excel("data-raw/meta.xlsx", sheet = "ratios"))
row.names(ratio) <- ratio$Numerator_Code

for (i in variables_stock){
  a <- ratio[i, "Denominator_Code"]
  b <- ratio[i, "Coeff"]
  c <- sum(df2011[, i], na.rm = TRUE)
  d <- sum(df2011[, a], na.rm = TRUE)
  moyfr <- round(b * c/d, 2)
  e <- sum(pourmoyennesCAH[, i], na.rm = TRUE)
  f <- sum(pourmoyennesCAH[, a], na.rm = TRUE)
  moyCFus <- round(b * e/f, 2)
  g <- ratio[i, "CODE"]
  h <- c(g, moyfr, moyCFus)
  moyennesCAH <- rbind(moyennesCAH, h, stringsAsFactors= FALSE)
  rm( a, b, c, d, e, f, g, h, moyfr, moyCFus)
}

colnames(moyennesCAH) <- c("Variable", "France", "CommunesFusionnantes")
moyennesCAH$Variable <- as.factor(moyennesCAH$Variable)
moyennesCAH$France <- as.numeric(moyennesCAH$France)
moyennesCAH$CommunesFusionnantes <- as.numeric(moyennesCAH$CommunesFusionnantes)
moyennesCAH$DiffComFusComFr <- moyennesCAH$CommunesFusionnantes - moyennesCAH$France

aa<- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = France), stat = "identity") + coord_flip()
bb <- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = CommunesFusionnantes), stat = "identity") + coord_flip()
cc <- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = DiffComFusComFr), stat = "identity") + coord_flip() + ylab(paste0("Différence Communes fusionnantes - ", "France"))
cowplot::plot_grid(aa, bb, cc, ncol = 1, nrow = 3)

# Pour variables politiques
# Calculs des valeurs moyennes (pour comparaison)
# On identifie les variables en stock
VarCAHBrutes <- stringr::str_replace(variables_elect, "prct_insc", "nbre_voix")
# Quelques cas particuliers
VarCAHBrutes <- stringr::str_replace(VarCAHBrutes, "Abst_nbre_voix", "Abst_nbre")
VarCAHBrutes <- stringr::str_replace(VarCAHBrutes, "BlcsNuls_nbre_voix", "BlcsNuls_nbre")

elect <- substr(VarCAHBrutes[1], start = 1, stop = 9)

# Sélection données utiles
pourmoyennesCAH <- na.omit(subset(df2011, COM_NOUV == "OUI"))
somme_inscr_CFus <- sum(pourmoyennesCAH[, paste0(elect, "_Inscr_nbre")], na.rm = TRUE)
somme_inscr_ttes_com <- sum(na.omit(df2011[, paste0(elect, "_Inscr_nbre")]))

moyennesCAH <- data.frame()

i <- VarCAHBrutes[3]

for (i in VarCAHBrutes){
  somme_ttes_com <- sum(df2011[, i], na.rm = TRUE)
  moy_ttes_com <- round(100*somme_ttes_com/somme_inscr_ttes_com, 2)
  somme_CFus <- sum(pourmoyennesCAH[, i], na.rm = TRUE)
  moy_Cfus <- round(100*somme_CFus/somme_inscr_CFus, 2)
  variable <- variables_elect[ which(VarCAHBrutes == i) ]
  ligne <- c(variable, moy_ttes_com, moy_Cfus)
  moyennesCAH <- rbind(moyennesCAH, ligne, stringsAsFactors= FALSE)
  
}
rm (somme_inscr_CFus, somme_inscr_ttes_com, somme_ttes_com, moy_ttes_com, somme_CFus, moy_Cfus, ligne, variable)

colnames(moyennesCAH) <- c("Variable", "Ensemble_étudié", "CommunesFusionnantes")
moyennesCAH$Variable <- as.factor(moyennesCAH$Variable)
moyennesCAH$Ensemble_étudié <- as.numeric(moyennesCAH$Ensemble_étudié)
moyennesCAH$CommunesFusionnantes <- as.numeric(moyennesCAH$CommunesFusionnantes)
moyennesCAH$DiffComFusComFr <- moyennesCAH$CommunesFusionnantes - moyennesCAH$Ensemble_étudié

aa<- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = Ensemble_étudié), stat = "identity") + coord_flip() + ylab("France")
bb <- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = CommunesFusionnantes), stat = "identity") + coord_flip() + ylab("Communes fusionnantes")
cc <- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = DiffComFusComFr), stat = "identity") + coord_flip() + ylab(paste0("Différence Communes fusionnantes - ", "France"))
cowplot::plot_grid(aa, bb, cc, ncol = 1, nrow = 3)

7 Régression logistique classique

Cf. par exemple https://larmarange.github.io/analyse-R/regression-logistique.html.

On fait le choix d’analyser ces données sous l’angle d’une régression logistique. Il s’agit d’observer si des variables augmentent ou diminuent la probabilité d’avoir une situation (fusion ou non). (Les tests de régression logistique sont utilisées parfois pour observer si des caractéristiques, qui seraient donc alors des comorbidités, augmentent ou non la probabilité d’avoir un décès.)

7.1 Variables explicatives envisagées

Les variables explicatives envisagées sont celles qui ont été pointées comme caractérisant davantage les communes nouvelles lors des études univariées, bivariées, ou lors de catégorisations Bideau and Ysebaert (2022).

Le choix des variables explicatives envisagées est déterminant, comme cela a pu être montré par exemple par Afsa (2016). En effet, le modèle logit présente la probabilité d’estimer la probabilité qu’un évènement se produise (dans notre cas, fusion ou non) mais il seulement en prenant en compte les données présentes dans le modèle. Dans le cas présenté par Afsa (2016), si on n’intègre pas les CSP des parents, le fait d’être en ZEP a un impact très négatif sur l’orientation en seconde. Si on intègre ces CSP, la significativité de la ZEP disparaît presque intégralement, si on ajoute le niveau en 6e, le fait d’être en ZEP redevient significatif mais un facteur favorable pour cette orientation.

# Si l'on souhaite avoir des noms de variables plus explicites, il faut ajouter des étiquettes des variables avec var_label de l'extension labelled. Par exemple :
# var_label(d$sport) <- "Pratique du sport ?"

7.2 Modalités de la variable d’intérêt

Il est d’abord pertinent de vérifier que ce qui va être perçu, dans notre variable d’intérêt (ici la fusion ou non) comme la modalité de référence est bien la situation normale, la plus fréquente. Dans notre cas, la très grandes majorité des communes françaises n’a pas connu de fusion, c’est donc la non-fusion qui est notre modalité de référence.

# Nécessité de passer la variable d'intérêt en facteur pour la fonction glm suivante
df_reg$COM_NOUV <- as.factor(df_reg$COM_NOUV)
class(df_reg$COM_NOUV)

## [1] "factor"

freq(df_reg$COM_NOUV)

##         n    % val%
## NON 33537 92.9 92.9
## OUI  2571  7.1  7.1

# Vérification que la variable la plus fréquente est bien la première (doit être la variable de référence)
freq(df_reg$COM_NOUV, sort = "dec")

##         n    % val%
## NON 33537 92.9 92.9
## OUI  2571  7.1  7.1

row.names(freq(df_reg$COM_NOUV, sort = "dec"))

## [1] "NON" "OUI"

ordre <- row.names(freq(df_reg$COM_NOUV, sort = "dec"))
df_reg$COM_NOUV <- factor (df_reg$COM_NOUV, ordre)

# Si on veut forcer dans un sens :
# df_reg$COM_NOUV <- relevel (df_reg$COM_NOUV, "NON", "OUI")

df_reg$ChefLieu <- as.factor(df_reg$ChefLieu)
freq(df_reg$ChefLieu, sort = "dec")

##       n    % val%
## O 34339 95.1 95.1
## N  1769  4.9  4.9

row.names(freq(df_reg$ChefLieu, sort = "dec"))

## [1] "O" "N"

ordre <- row.names(freq(df_reg$ChefLieu, sort = "dec"))
df_reg$ChefLieu <- factor (df_reg$ChefLieu, ordre)

7.3 Mise en œuvre de la régression logistique

Choix des variables qui doivent éviter celles pointées comme auto-corrélées plus haut. On affiche également le VIF pour limiter au maximum les biais d’interprétation.

L’approche la plus classique consiste à examiner les facteurs d’inflation de la variance (FIV) ou variance inflation factor (VIF) en anglais. Les FIV estimenent de combien la variance d’un coefficient est augmentée en raison d’une relation linéaire avec d’autres prédicteurs. Ainsi, un FIV de 1,8 nous dit que la variance de ce coefficient particulier est supérieure de 80 % à la variance que l’on aurait dû observer si ce facteur n’est absolument pas corrélé aux autres prédicteurs.

Si tous les FIV sont égaux à 1, il n’existe pas de multicolinéarité, mais si certains FIV sont supérieurs à 1, les prédicteurs sont corrélés. Il n’y a pas de consensus sur la valeur au-delà de laquelle on doit considérer qu’il y a multicolinéarité. Certains auteurs, comme Paul Allison, disent regarder plus en détail les variables avec un FIV supérieur à 2,5. D’autres ne s’inquiètent qu’à partir de 5. Il n’existe pas de test statistique qui permettrait de dire s’il y a colinéarité ou non. (https://larmarange.github.io/analyse-R/multicolinearite.html)

# Si on veut importer directement les données du modèle réalisé précédemment
# load("data/refdata_modele.Rdata")




# Les variables explicatives sont mentionnées juste après la variable d'intérêt (ici, "COM_NOUV")
reg <- glm(COM_NOUV ~ 
# reg <- glm(ChefLieu ~ 
             # dgf_prct + perso_prct + equip_prct + dette_prct,
             # dgf_prct + perso_prct + equip_prct,
             # P09_RETR1564_RT + C09_ACT1564_Agr_RT + P09_POP3044Y_RT + C09_ACTOCC_OUT_RT + P09_CHOM1564_RT + dgf_prct + perso_prct + charge_prct,
             # P09_RETR1564_RT + C09_ACT1564_Agr_RT + P09_POP3044Y_RT + C09_ACTOCC_OUT_RT + P09_CHOM1564_RT + dgf_prct + perso_prct,
             # P09_RETR1564_RT + C09_ACT1564_Agr_RT + C09_ACTOCC_OUT_RT + P09_CHOM1564_RT,
             # P09_ETUD1564_RT + C09_ACT1564_Ouvr_RT + C09_ACTOCC_OUT_RT + P09_CHOM1564_RT,
             # C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT + dgf_prct,
             # C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali + dgf_prct_quali, # Effets faibles mais significatifs
             # P09_ETUD1564_RT_quali + C09_ACT1564_Ouvr_RT_quali + C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali,
             # C09_ACTOCC_OUT_RT_quali,
             # C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali + P09_POP1529Y_RT_quali + P09_POP_quali + PR2012_T1_Abst_prct_insc_quali + PR2012_T1_prct_insc_LE.PEN_quali + CATAEU2010 + CODE_DEPT,
             # P09_ETUD1564_RT_quali + C09_ACT1564_Ouvr_RT_quali + C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali + P09_POP1529Y_RT_quali + P09_POP_quali + PR2012_T1_Abst_prct_insc_quali + PR2012_T1_prct_insc_LE.PEN_quali + dgf_prct_quali + perso_prct_quali + CATAEU2010,
             # P09_ETUD1564_RT_quali + C09_ACT1564_Ouvr_RT_quali + C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali + P09_POP1529Y_RT_quali + P09_POP_quali + PR2012_T1_Abst_prct_insc_quali + PR2012_T1_prct_insc_LE.PEN_quali + dgf_prct_quali + perso_prct_quali + CATAEU2010 + superficie_quali,
             # P09_ETUD1564_RT_quali + C09_ACT1564_Ouvr_RT_quali + C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali + P09_POP1529Y_RT_quali + P09_POP_quali + PR2012_T1_Abst_prct_insc_quali + PR2012_T1_prct_insc_LE.PEN_quali + dgf_prct_quali + perso_prct_quali + CATAEU2010 + superficie_quali + CODE_DEPT,
              # CATAEU2010 + P09_POP_quali,
              # CATAEU2010 + RPI,
              # CATAEU2010 + `École élémentaire_quantile`,
              # CATAEU2010 + `École élémentaire_quantile` + `Médecin omnipraticien_quantile`,
              # CATAEU2010 + `Médecin omnipraticien_quantile`,

# Variables ressortant des pas à pas, classées de la plus significative à la moins significative
             # CODE_DEPT + CATAEU2010 + P11_IMP_NET_RT_quali + C09_EMP_CONC_RT_quali + PR2012_T1_prct_insc_LE.PEN_quali  + PotFin4taxes_prct_quali  + perso_prct_quali  + res1_prct_quali + Supermarché_quantile + P09_POP_quali + `Médecin omnipraticien_quantile` + P09_POP1529Y_RT_quali + P11_POT_FIN_RT_quali + P09_POP0014Y_RT_quali + dette_prct_quali + P11_FoyFisc_Imp_RT_quali + PR2012_T1_prct_insc_BAYROU_quali + P09_POP4559Y_RT_quali + PR2012_T1_Abst_prct_insc_quali,

# Variables s'inspirant des pas à pas (test1)
             # CODE_DEPT + P09_POP_quali + CATAEU2010 + `Médecin omnipraticien_quantile` + PotFin4taxes_prct_quali + Supermarché_quantile + C09_EMP_CONC_RT_quali + PR2012_T1_prct_insc_LE.PEN_quali + P09_POP1529Y_RT_quali + dgf_prct_quali,
# Variables s'inspirant des pas à pas, sans départements (test2)
             P09_POP_quali + CATAEU2010 + `Médecin omnipraticien_quantile` + PotFin4taxes_prct_quali + Supermarché_quantile + C09_EMP_CONC_RT_quali + PR2012_T1_prct_insc_LE.PEN_quali + P09_POP1529Y_RT_quali + dgf_prct_quali,

             # P09_POP_quali + CATAEU2010 + `Médecin omnipraticien_quantile` + PotFin4taxes_prct_quali + Supermarché_quantile + C09_EMP_CONC_RT_quali + PR2012_T1_prct_insc_LE.PEN_quali + P09_POP1529Y_RT_quali + dgf_prct_quali + PR2012_T1_Abst_prct_insc_quali,

# Quelques variables simples
             # P09_POP_quali + superficie_quali + CATAEU2010 + CODE_DEPT,
             # P09_POP_quali + superficie_quali + CATAEU2010,
             # P09_POP_quali +  CATAEU2010 + CODE_DEPT,

           data = df_reg, family = binomial(logit), na.action = na.exclude)

options(max.print=700)
summary(reg)
# exp(coef(reg))
# exp(confint(reg))
# exp(cbind(coef(reg), confint(reg)))
# odds.ratio(reg)

# Pour vérifier les VIF
car::vif(reg)
which(car::vif(reg)[,3] > 1.9)
# test <- as.data.frame(car::vif(reg))

# Présentation d'un tableau plus propre
# tbl_regression(reg, exponentiate = TRUE)
add_vif(tbl_regression(reg, exponentiate = TRUE)) # Avec affichage des VIF

# Représentations graphiques du modèle
# ggcoef_model(reg, exponentiate = FALSE)
forest_model(reg)

# Représentation graphique des effets
# plot(allEffects(reg))

# plot_grid(plotlist = plot(ggeffect(reg)))
# Pour n'afficher qu'une variable :
# plot(ggeffect(reg, "_nom_de_la_variable"))

screenreg(reg, stars = c(0.01, 0.05, 0.1))

# texreg(reg, stars = c(0.01, 0.05, 0.1)) # Pour export LaTeX

7.4 Lecture du modèle

Les odds ratio (littéralement rapport des cotes, ou rapport des chances ou rapport des risques relatifs) permettent de décrire la significativité pratique. Celle-ci désigne l’impact d’une variable. Il ne faut pas la confondre avec la significativité statistique, qui décrit « le degré de certitude avec lequel on peut affirmer qu’une variable influe » (Afsa (2016) p. 39).

Il développe : « L’odds ratio attaché à la variable zep est égal à 1,631, avec [1.416 , 1.878] comme intervalle de confiance à 95%. On l’a vu en section I.6.b, cela signifie précisément que la chance relative de passer en seconde générale est environ 1,6 fois plus élevé pour un enfant en ZEP que pour un enfant hors ZEP, conditionnellement aux facteurs pris en compte dans le modèle (i.e. à sexe, âge, milieu social et niveau en sixième fixés). Il est entendu que la chance relative est un rapport de probabilités : c’est la probabilité de passer en seconde générale rapportée à celle de ne pas y passer. L’odds ratio n’est donc pas un rapport de deux probabilités, mais un rapport de rapports de probabilités. Il ne faut surtout pas dire que le fait d’être en ZEP multiplie par 1,6 la probabilité de passer en seconde générale, à mêmes caractéristiques observées. On verra plus loin que ce résultat est, en ces termes, complètement faux. » (p. 71)

À noter que, par exemple pour Afsa (2016), l’expression « toutes choses égales par ailleurs » est à éviter. En effet, ce n’est le cas que si les variables étudiées permettent de décrire parfaitement le phénomène. « les résultats des estimations sont conditionnels à la liste des variables x introduites dans le modèle, c’est-à-dire qu’ils dépendent des variables introduites » (Afsa (2016) p. 53). Si des variables, non mesurées, influent à la fois sur les variables du modèle et la variable étudiée, cela brouille l’analyse. Ainsi, pour les analyses présentées ici, c’est plutôt l’expression “Toutes choses égales quant aux autres variables introduites” qui est à utiliser.

Un Odd ratio égal à 1 signifie que la variable d’intérêt (la fusion ici) est indépendante de la variable explicative envisagée. Si l’OR est supérieur à 1, une commune ayant une variable explicative plus élevée ont plus de chance de fusionner, s’il est inférieur à 1, les communes avec une variable explicative plus faible ont moins de chance de fusionner.

7.5 Sauvegarde données modèle

# try(save(df_reg, reg, file = "data/modeles/refdata_modele_test_ttes_communes.Rdata"))
# try(save(df_reg, reg, file = "data/modeles/refdata_modele_test2.Rdata"))

load("data/modeles/refdata_modele_test2.Rdata")


# try(save(test, file = "data/modeles/refdata_modele_test1_tableau_resultats.Rdata"))
# load("data/modeles/refdata_modele_test1_tableau_resultats.Rdata")

7.6 Cartographie de l’effet des départements

Pour les modèles intégrant les appartenances départementales, il n’est pas aisé de représenter les résultats dans un tableau. Choix, donc, de réaliser une carte.

# test <- as.data.frame(tbl_regression(reg, exponentiate = TRUE))

# try(save(test, file = "data/modeles/refdata_modele_test1_tableau_resultats.Rdata"))
load("data/modeles/refdata_modele_test1_tableau_resultats.Rdata")

tableau_departements <- test [3:92,c(1, 2, 4)]
colnames(tableau_departements) <- c("CODE_DEPT", "OR", "pvalue")
tableau_departements$pvalue[tableau_departements$pvalue == "<0.001"] <- 0.0001
tableau_departements$pvalue[tableau_departements$pvalue == ">0.9"] <- 1
# On passe les variables en numérique
tableau_departements <- as.data.frame(apply(tableau_departements, 2, as.numeric))

# Choix d'un seuil à partir duquel on ne représente plus les départements
seuil_choisi <- 0.1
tableau_departements <- subset(tableau_departements, tableau_departements$pvalue <= seuil_choisi)

tableau_departements <- merge(dep, tableau_departements, by = "CODE_DEPT", all.x = TRUE)
par(mar=c(0,2,1.2,2))

layoutLayer(title = " ",
            author = "G. Bideau, 2024",
            sources = paste0( "Les départements non renseignés ont\nune p-value supérieure à ", seuil_choisi, ".\n\n",
                              "Sources : INSEE, IGN, 2024"), extent = dep,
            frame = FALSE,
            col = "white",
            tabtitle = FALSE)

plot(st_geometry(dep), col = NA, add = TRUE)

propSymbolsChoroLayer(x = tableau_departements, var = "OR", var2 = "pvalue",
                      col = rev(carto.pal(pal1 = "red.pal", n1 = 4)), # On inverse l'ordre pour avoir du rouge près de zéro
                      breaks = c(0, 0.001, 0.01, 0.05, 0.1),
                      inches = 0.2, method = "q4",
                      border = "grey50", lwd = 1,
                      legend.var.pos = "topright", 
                      legend.var2.pos = "topleft",
                      legend.var2.values.rnd = 4,
                      legend.var2.title.txt = "Valeur de la p-value",
                      legend.var.title.txt = "Odd Ratio",
                      legend.var.style = "e",
                      legend.title.cex = 0.9,
                      legend.values.cex = 0.8,
                      add = TRUE)

7.7 Identifier les variables ayant un effet significatif

Cf. https://larmarange.github.io/analyse-R/regression-logistique.html#identifier-les-variables-ayant-un-effet-significatif : ” Pour tester l’effet global sur un modèle, on peut avoir recours à la fonction drop1. Cette dernière va tour à tour supprimer chaque variable du modèle et réaliser une analyse de variance (ANOVA, voir fonction anova) pour voir si la variance change significativement.”

drop1(reg, test = "Chisq")

## Single term deletions
## 
## Model:
## COM_NOUV ~ CODE_DEPT + P09_POP_quali + CATAEU2010 + `Médecin omnipraticien_quantile` + 
##     PotFin4taxes_prct_quali + Supermarché_quantile + C09_EMP_CONC_RT_quali + 
##     PR2012_T1_prct_insc_LE.PEN_quali + P09_POP1529Y_RT_quali + 
##     dgf_prct_quali
##                                  Df Deviance   AIC     LRT  Pr(>Chi)    
## <none>                                 14038 14300                      
## CODE_DEPT                        92    17002 17080 2964.75 < 2.2e-16 ***
## P09_POP_quali                     4    14089 14343   51.41 1.832e-10 ***
## CATAEU2010                        8    14110 14356   72.80 1.359e-12 ***
## `Médecin omnipraticien_quantile`  4    14068 14322   30.00 4.899e-06 ***
## PotFin4taxes_prct_quali           4    14058 14312   20.43 0.0004100 ***
## Supermarché_quantile              2    14055 14313   16.98 0.0002053 ***
## C09_EMP_CONC_RT_quali             4    14053 14307   15.66 0.0035126 ** 
## PR2012_T1_prct_insc_LE.PEN_quali  4    14052 14306   14.09 0.0070103 ** 
## P09_POP1529Y_RT_quali             4    14053 14307   15.55 0.0036846 ** 
## dgf_prct_quali                    4    14040 14294    2.28 0.6847227    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kable(as.data.frame(drop1(reg, test = "Chisq"))) # Pour identifier les variables ayant un effet significatif

	Df	Deviance	AIC	LRT	Pr(>Chi)
	NA	14037.59	14299.59	NA	NA
CODE_DEPT	92	17002.34	17080.34	2964.747857	0.0000000
P09_POP_quali	4	14089.00	14343.00	51.410589	0.0000000
CATAEU2010	8	14110.39	14356.39	72.799731	0.0000000
`Médecin omnipraticien_quantile`	4	14067.59	14321.59	29.997805	0.0000049
PotFin4taxes_prct_quali	4	14058.03	14312.03	20.433467	0.0004100
Supermarché_quantile	2	14054.58	14312.58	16.982134	0.0002053
C09_EMP_CONC_RT_quali	4	14053.25	14307.25	15.659018	0.0035126
PR2012_T1_prct_insc_LE.PEN_quali	4	14051.68	14305.68	14.090965	0.0070103
P09_POP1529Y_RT_quali	4	14053.14	14307.14	15.551106	0.0036846
dgf_prct_quali	4	14039.87	14293.87	2.278299	0.6847227

test_anova <- Anova(reg)
kable(as.data.frame(test_anova))

	LR Chisq	Df	Pr(>Chisq)
CODE_DEPT	2964.747857	92	0.0000000
P09_POP_quali	51.410589	4	0.0000000
CATAEU2010	72.799731	8	0.0000000
`Médecin omnipraticien_quantile`	29.997805	4	0.0000049
PotFin4taxes_prct_quali	20.433467	4	0.0004100
Supermarché_quantile	16.982134	2	0.0002053
C09_EMP_CONC_RT_quali	15.659018	4	0.0035126
PR2012_T1_prct_insc_LE.PEN_quali	14.090965	4	0.0070103
P09_POP1529Y_RT_quali	15.551106	4	0.0036846
dgf_prct_quali	2.278299	4	0.6847227

Anova(reg, type = "II")

## Analysis of Deviance Table (Type II tests)
## 
## Response: COM_NOUV
##                                  LR Chisq Df Pr(>Chisq)    
## CODE_DEPT                         2964.75 92  < 2.2e-16 ***
## P09_POP_quali                       51.41  4  1.832e-10 ***
## CATAEU2010                          72.80  8  1.359e-12 ***
## `Médecin omnipraticien_quantile`    30.00  4  4.899e-06 ***
## PotFin4taxes_prct_quali             20.43  4  0.0004100 ***
## Supermarché_quantile                16.98  2  0.0002053 ***
## C09_EMP_CONC_RT_quali               15.66  4  0.0035126 ** 
## PR2012_T1_prct_insc_LE.PEN_quali    14.09  4  0.0070103 ** 
## P09_POP1529Y_RT_quali               15.55  4  0.0036846 ** 
## dgf_prct_quali                       2.28  4  0.6847227    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

reg %>%
  tbl_regression(exponentiate = TRUE) %>%
  add_global_p(type = "II")

Caractéristique	OR¹	95% IC¹	p-valeur
CODE_DEPT			<0,001
62	—	—
02	2,93	1,39 – 6,70
80	1,87	0,83 – 4,45
76	6,51	3,32 – 14,3
57	1,36	0,54 – 3,45
14	38,5	20,7 – 82,0
21	2,14	0,95 – 5,11
60	3,17	1,48 – 7,35
27	20,0	10,6 – 42,9
59	0,72	0,19 – 2,23
51	2,31	1,02 – 5,56
50	42,0	22,4 – 89,7
25	6,52	3,27 – 14,5
54	0,86	0,26 – 2,52
31	1,07	0,36 – 3,00
71	2,08	0,89 – 5,09
24	13,5	6,99 – 29,3
64	0,66	0,18 – 2,04
70	1,74	0,69 – 4,43
39	11,5	5,94 – 25,2
33	2,46	1,05 – 6,02
38	8,12	4,05 – 18,1
67	3,71	1,72 – 8,67
88	2,42	1,01 – 5,98
77	2,31	0,94 – 5,79
61	32,4	17,1 – 69,5
55	0,20	0,01 – 1,08
65	1,36	0,48 – 3,70
17	3,06	1,36 – 7,33
63	2,27	0,95 – 5,63
08	5,15	2,45 – 11,8
32	0,74	0,20 – 2,30
89	8,98	4,51 – 19,9
11	1,95	0,75 – 5,05
52	4,60	2,10 – 10,8
10	0,67	0,15 – 2,28
01	11,7	5,89 – 26,1
16	15,8	8,11 – 34,6
28	12,2	6,19 – 26,8
68	5,37	2,46 – 12,6
72	9,35	4,61 – 21,0
22	11,0	5,46 – 24,7
26	2,72	1,11 – 6,84
49	168	88,0 – 362
30	1,26	0,34 – 3,92
35	11,5	5,68 – 26,0
34	0,62	0,09 – 2,42
46	10,2	5,04 – 22,9
07	1,93	0,68 – 5,26
45	3,73	1,56 – 9,26
09	2,36	0,91 – 6,15
40	1,59	0,53 – 4,49
42	2,27	0,80 – 6,15
81	3,79	1,61 – 9,31
03	1,37	0,42 – 4,00
47	0,00	0,00 – 0,00
58	0,55	0,08 – 2,16
12	7,78	3,70 – 17,9
73	17,5	8,80 – 38,7
79	25,6	13,1 – 56,3
74	10,0	4,69 – 23,2
41	12,0	5,85 – 27,1
69	18,1	8,88 – 40,9
18	1,23	0,33 – 3,82
19	4,32	1,85 – 10,6
85	17,6	8,80 – 39,4
29	4,30	1,77 – 10,8
86	9,02	4,25 – 20,9
37	3,33	1,28 – 8,66
53	12,9	6,32 – 29,2
78	2,65	0,88 – 7,50
15	6,81	3,11 – 16,1
23	2,59	0,95 – 6,91
43	1,33	0,36 – 4,13
56	8,15	3,76 – 19,1
36	3,67	1,45 – 9,40
66	0,00	0,00 – 0,00
44	13,8	6,36 – 32,3
04	1,65	0,44 – 5,16
87	4,30	1,64 – 11,2
91	2,70	0,72 – 8,45
82	0,00	0,00 – 0,00
48	29,0	14,4 – 65,1
95	2,57	0,69 – 8,03
05	8,81	3,94 – 21,1
06	0,00	0,00 – 0,00
83	0,00	0,00 – 0,00
84	0,00	0,00 – 0,00
13	0,00	0,00 – 0,00
90	1,96	0,29 – 7,76
94	0,00	0,00 – 0,00
93	0,00	0,00 – 0,00
92	0,00	0,00 – 0,00
P09_POP_quali			<0,001
Q3	—	—
Q1	1,05	0,90 – 1,23
Q2	1,03	0,90 – 1,18
Q4	0,73	0,63 – 0,84
Q5	0,44	0,35 – 0,56
CATAEU2010			<0,001
112	—	—
400	1,36	1,17 – 1,58
300	1,18	1,03 – 1,35
120	1,16	0,99 – 1,36
111	0,44	0,33 – 0,57
221	1,28	0,96 – 1,70
212	1,15	0,86 – 1,52
222	0,88	0,56 – 1,32
211	0,85	0,53 – 1,32
Médecin omnipraticien_quantile			<0,001
[0,1)	—	—
[1,2)	1,26	1,05 – 1,53
[2,4)	1,72	1,39 – 2,12
[4,7)	1,52	1,12 – 2,05
[7,542)	1,02	0,68 – 1,53
PotFin4taxes_prct_quali			<0,001
Q3	—	—
Q1	0,76	0,65 – 0,90
Q2	0,79	0,68 – 0,92
Q4	0,87	0,76 – 1,00
Q5	1,02	0,89 – 1,18
Supermarché_quantile			<0,001
[0,1)	—	—
[1,2)	1,31	1,05 – 1,62
[2,55)	1,89	1,39 – 2,55
C09_EMP_CONC_RT_quali			0,004
Q3	—	—
Q1	1,22	1,05 – 1,43
Q2	1,00	0,87 – 1,16
Q4	1,01	0,87 – 1,17
Q5	1,21	1,03 – 1,42
PR2012_T1_prct_insc_LE.PEN_quali			0,007
Q3	—	—
Q1	1,02	0,88 – 1,18
Q2	1,03	0,90 – 1,17
Q4	1,04	0,90 – 1,20
Q5	0,78	0,66 – 0,92
P09_POP1529Y_RT_quali			0,004
Q3	—	—
Q1	0,79	0,68 – 0,92
Q2	0,81	0,70 – 0,93
Q4	0,87	0,76 – 1,00
Q5	0,98	0,85 – 1,13
dgf_prct_quali			0,7
Q3	—	—
Q1	0,96	0,82 – 1,13
Q2	1,01	0,87 – 1,17
Q4	0,92	0,80 – 1,06
Q5	0,92	0,79 – 1,07
¹ OR = rapport de cotes, IC = intervalle de confiance

On observe que les variables des navetteurs et de l’abstention sont potentiellement peu importantes pour le modèle.

8 Étude des résidus (y compris avec indice de Moran pour une auto-corrélation spatiale des résidus ?)

Non-stationnarité spatiale : lorsque les paramètres permettant d’expliquer un phénomènes sont différents en fonction des espaces. Bonne figure p. 185 (Feuillet et al, 2019) où en regardant deux groupes distincts on aperçoit un type de relation alors que lorsqu’on regarde les deux en même temps, on imagine une autre relation statistique. C’est le paradoxe de Simpson (qui n’est pas forcément situé).

Il est donc intéressant voire important de cartographier ou de représenter graphiquement les résidus. Cela permet potentiellement de valider ou invalider le modèle mais également de réfléchir aux espaces où ce modèle fonctionne le mieux ou le moins bien.

8.1 Étude résidus régression linéaire (peu pertinent dans le cas d’un logit)

Cf. R et Espace section 5.2.3 et particulièrement p. 91 sq pour aperçu des résidus : « La fonction plot() appliquée aux résultats d’un modèle de régression linéaire obtenus par la fonction lm() permet de représenter les quatre principales hypothèses au cœur de ce modèle :

— la normalité des résidus par rapport aux valeurs prédites (en haut à gauche de l’image)

— la normalité globale des résidus (en haut à droite de l’image)

— la corrélation entre les valeurs de la variable explicative et le carré des résidus standardisés (en bas à gauche de l’image)

— l’existence de valeurs extrêmes altérant l’estimation des paramètres (en bas à droite de l’image) »

Autre source très utilisée ci-dessous : https://sites.google.com/site/rgraphiques/4--stat/machine-learning-biostatistiques-analyse-de-donn%C3%A9es/r%C3%A9gressions-lin%C3%A9aires-avec-r?pli=1

# On veut éviter les phénomènes d'auto-corrélation : les résidus augmenteraient ensemble dans une zone donnée et seraient liés les uns aux résidus précédents ou suivants (cas typiques lorsque X est un enregistrement du temps). Vérification par le test de Durbin-Watson (indépendance : p-value > 0,05)
durbinWatsonTest(reg)
dwtest(reg)

# Visualisation des résidus
par(mfrow = c(2,2))
plot(reg)
par(mfrow = c(1,1))

# Les résidus doivent être répartis approximativement selon une courbe de Gauss
hist(residuals(reg),col="yellow",freq=F)

# Il faut, pour que la régression soit pertinente, observer un alignement entre les quantiles des résidus et les quantiles théoriques, donc globalement une ligne droite
plot(reg, which = 2)
# On peut vérifier avec le test de normalité de Shapiro-Wilk (normalité si p-value supérieure à 0,05)
shapiro.test(residuals(reg))


# Vérification que la variance des résidus est constante (distribution homogène)
# Pour que le modèle soit pertinent, il faut une courbe rouge plane et que les valeurs ne se regroupent pas.
plot(reg, which = 3)
# On peut vérifier avec le test de Breush-Pagan (homogénéité si p-value > 0,05)
# ncvTest(reg) # seulement pour objet lm
bptest(reg)
# Ou test de Goldfeld-Quandt (homogénéité si p-value > 0,05)
gqtest(reg)

plot(reg, which = 5)

8.2 Étude résidus logit

Cf. https://pmarchand1.github.io/ECL7102/notes_cours/9-Regression_logistique.pdf#page=9

library(arm)

binnedplot(fitted(reg), residuals(reg, type = "response"))

# Graphique prédiction linéaire/résidus
# Cf. https://lrouviere.github.io/doc_cours/poly_logistique.pdf#page=49
res_dev <- residuals(reg) #residus de deviance
res_pear <- residuals(reg,type="pearson") #residus de Pearson
res_dev_stand <- rstandard(reg) #residu de deviance standardises 
H <- influence(reg)$hat #diagonale de la hat matrix
res_pear_stand <- res_pear/sqrt(1-H) #residu de Pearson standardises
plot(rstudent(reg),type="p",cex=0.5,ylab="Résidus studentisés par VC")
abline(h=c(-2,2))

# Résidus partiels (résidus pour une variable)
# Cf. https://lrouviere.github.io/doc_cours/poly_logistique.pdf#page=50
residpartiels<-resid(reg,type="partial")
prov<-loess(residpartiels[,"P09_POP_quali"]~df_reg$P09_POP_quali)
ordre<-order(df_reg$P09_POP_quali)
plot(df_reg$P09_POP_quali,residpartiels[,"P09_POP_quali"],type="p",cex=0.5,xlab="",ylab="")
matlines(df_reg$P09_POP_quali[ordre],predict(prov)[ordre])
abline(lsfit(df_reg$P09_POP_quali,residpartiels[,"P09_POP_quali"]),lty=2)

8.3 Cartographie des résidus

# residuals (reg)
# On intègre les résidus aux données
df_reg$residus <- residuals (reg)
par(mar=c(8,0,8,0), mfrow = c(1,2))


choroLayer(x = df_reg , var = "residus", 
           method = "quantile", nclass = 6,
           col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           border = NA,
           colNA = "white",
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Résidus",
           # add = TRUE
           )

plot(st_geometry(dep), col = NA,
     add = TRUE
     )

layoutLayer(title = "Cartographie des résidus du modèle",
            author = "Auteur : G. Bideau, 2024.",
            sources = "Sources : INSEE, 2024.", col = "white", coltitle = "black")

typoLayer(x = df_reg, var = "COM_NOUV",
          legend.values.order = c("OUI", "NON"),
          legend.pos = "topleft",
          col=c("red","blue"),
          legend.title.txt = "Communes fusionnantes",
          border = NA)

layoutLayer(title = "Participation à une commune nouvelle 2012-2024(01)",
            author = "Auteur : G. Bideau, 2024.",
            sources = "Sources : INSEE, 2024.", col = "white", coltitle = "black")

par(mar = c(0,0,0,0), mfrow = c(1,1))

knitr::opts_chunk$set(echo=FALSE, # Afficher ou non le code R dans le document
                      eval=FALSE, # Exécuter ou non le code R à la compilation
                      include=TRUE, #   Inclure ou non le code R et ses résultats dans le document
                      # results “hide”/“asis”/“markup”/“hold”   Type de résultats renvoyés par le bloc de code
                      warning=TRUE, # Afficher ou non les avertissements générés par le bloc
                      message=TRUE, # Afficher ou non les messages générés par le bloc
                      cache=TRUE) # Utiliser le cache pour accélerer les knits.

8.4 Test de Moran sur auto-corrélation spatiale des résidus (package rgeoda)

9 Sélection step by step, Forward stepwise ou backward stepwise

Possibilité de sélectionner automatiquement les variables à choisir pour une régression logistique [feuillet2019] p. 174 :

régression pas à pas ascendante (on ajoute des variables à partir d’un tableau de données) : forward stepwise ;
régression pas à pas descendante (on enlève) : backward stepwise.

Guides pour mettre en place ces méthodes : https://www.statology.org/stepwise-regression-r/ pour une présentation un peu rapide ou (déroulé suivi ici) : https://rstudio-pubs-static.s3.amazonaws.com/205694_3b195f29e9504d23aeb483ff1ffafeba.html, en particulier à partir du point 3.4.

9.1 Sélection des données

9.2 Vérification des auto-corrélations et des manques de données

9.3 Régression pas à pas

9.3.1 Préparation données pour régression stepwise

9.3.2 Stepwise ascendante (forward)

9.3.3 Stepwise descendante (backward)

9.3.4 Stepwise dans les deux sens

Non joué pour l’instant car bloque ?

9.3.5 Stepwise F-test ascendante (forward)

9.3.6 Stepwise F-test descendante (backward)

9.3.7 Stepwise F-test deux sens

9.4 Tentative algorithme génétique

Section 3.5 de la page déjà citée : https://rstudio-pubs-static.s3.amazonaws.com/205694_3b195f29e9504d23aeb483ff1ffafeba.html#lalgorithme-genetique

Ne fonctionne pas (message d’erreur : “Erreur dans .subset2(x, i, exact = exact) : tentative de sélection de moins d’un élément dans get1index”).

Bibliographie

Afsa, Cédric. 2016. “Le Modèle Logit : Théorie Et Application - Documents de Travail - M2016/01.” INSEE. https://www.insee.fr/fr/statistiques/2022139.

Bideau, Gabriel. 2019. “Les communes nouvelles françaises (2010-2019) : Une réforme territoriale silencieuse.” Annales de Géographie 728 (4/2019): 57–85. https://www.cairn.info/revue-annales-de-geographie-2019-4-page-57.htm.

Bideau, Gabriel, and Ronan Ysebaert. 2022. “Les communes nouvelles françaises (2010-2020) : quels profils pour ces territoires du quotidien remaniés ?” L’Espace géographique 50 (2021/1-2): 45–66. https://doi.org/10.3917/eg.501.0045.

#CN_regressions Les communes nouvelles françaises (2012-2024) : proposition d’un modèle explicatif à l’aide de régressions logistiques

Gabriel Bideau

23 septembre 2024 13:15:03